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We recently obtained a quantum-Boltzmann-conserving classical dynamics by making a single change to the 
derivation of the ‘Classical Wigner’ approximation. Here, we show that the further approximation of this 
‘Matsubara dynamics’ gives rise to two popular heuristic methods for treating quantum Boltzmann time- 
correlation functions: centroid molecular dynamics (CMD) and ring-polymer molecular dynamics (RPMD). 
We show that CMD is a mean-field approximation to Matsubara dynamics, obtained by discarding (classi¬ 
cal) fluctuations around the centroid, and that RPMD is the result of discarding a term in the Matsubara 
Liouvillian which shifts the frequencies of these fluctuations. These findings are consistent with previous 
numerical results, and give explicit formulae for the terms that CMD and RPMD leave out. Copyright (2015) 
American Institute of Physics. This article may be downloaded for personal use only. Any other use requires 
prior permission of the author and the American Institute of Physics. The following article appeared in the 
Journal of Chemical Physics, 142, 191101 (2015) and may be found at http://dx.doi.org/10.1063/1.4921234 


I. INTRODUCTION 

Quantum Boltzmann time-correlation functions play a 
central role in chemical physics, and are (usually) impos¬ 
sible to calculate exactly. One promising approach is to 
treat the statistics quantally and the dynamics classically. 
The standard way to do this is to use the linearized semi- 
classical initial value representation (LSC-IVR or ‘classi¬ 
cal Wigner approximation’), 13 but this has the drawback 
of not satisfying detailed balance. Recently, 13 however, 
we found that a single change to the LSC-IVR derivation 
gives a classical dynamics which does satisfy detailed bal¬ 
ance. This modified version of LSC-IVR is called ‘Mat¬ 
subara dynamics’. 

We can summarise Matsubara dynamics as follows: At 
initial time, the quantum statistics gives rise to delocal¬ 
ized distributions in position which are smooth functions 
of imaginary time. If we constrain the LSC-IVR dynam¬ 
ics to conserve this smoothness (by including only the 
smooth ‘Matsubara’ modes—see Sec. II) we find that 
it satisfies detailed balance, and gives better ^reement 
than LSC-IVR with the exact quantum resulf.l^We sus¬ 
pect (but have not yet proved) that Matsubara dynam¬ 
ics reproduces the time-dependence of the exact Kubo- 
transformed time-correlation function up to order hP, and 
is thus the correct theory for describing quantum statis¬ 
tics and classical dynamics. 

Matsubara dynamics suffers from the sign problem and 
is hence impractical, but the findings just described sug¬ 
gest that it should be the starting point from which 
to make further approximations if one wishes to de¬ 
vise practical methods that combine quantum statistics 
with classical dynamics. Numerical tests in ref. (see 
also Fig. 1) showed that the popular centroid molec¬ 
ular dynamic^311 (CMD) and ring-polymer molecular 


dynamic^^H (RPMD) methods appear to be two such 
approximations. Here we confirm this, by deriving the 
terms that CMD and RPMD leave out from the Matsub¬ 
ara dynamics.^ 


II. SUMMARY OF MATSUBARA 
DYNAMICS 


Matsubara dynamics approximates the quantum 
Kubo-transformed time-correlation functiorP 
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and aM = hf^ ~ l)/2]!^- The position coordi¬ 

nates Q = {Q„}, with n = -{M - l)/2,..., (M - l)/2, 
are the M Matsubara modes, which describe closed paths 
q(r) that are smooth functions of the imaginary time t 
(= 0 — /h), where Qq is the centroid coordinate (see the 
Appendix); JdQ = J^^dQn, and P are similarly de¬ 
fined for momentum. The functions A(Q) and H(Q) are 
obtained from the operators A and B (see the Appendix), 
such that A = B = q gives A(Q) = B{Q) = Qo-^ The 
propagator 6"^“* contains the Matsubara Liouvillian 


(M-l)/2 

Cm = X 

n=-{M-l)/2 


Pu £ 

^ dQn 


dUniQ) 5 

dQn dPn 


Current address: Lab. fiir Physikalische Chemie, ETH Zurich, 
CH-8093 Zurich, Switzerlanc 
^^Corresponding author: scalO@cam.ac.uk 


1 


(4) 








in which the potential energy Um{Q,) is given in the Ap¬ 
pendix. The quantum Boltzmann distribution in Eq. ^ 
is complex, and contains the Matsubara Hamiltonian 

HMiP,Q) = ^+UM{Q) (5) 

2m 

and the phase 

(M-l)/2 

Om(P,Q)= ^ PntJnQ-n ( 6 ) 

n=-(M-l)/2 

where are the Matsubara frequencies = 2Tm/l3h. 
Matsubara dynamics is inherently classical (meaning 
that terms 0{h?) disappear from the quantum Liouvil- 
lian on decoupling the Matsubara modes, leaving Lm), 
and conserves the Hamiltonian Hm(P, Q) and the phase 
6m(P, Q), and thus satisfies detailed balance. 

Clearly Eq. (j^ suffers from the sign problem because of 
the phase 0 m(P,Q)- Let us make the coordinate trans¬ 
formation Pn = Pn — iniuinQ-n- This gives 
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is the ‘ring-polymer’ Hamiltonian familiar from quantum 
statistics.Q^®^ Equation Q is simply Eq. ^ in disguise, 
but at t = 0, we can use a standard contour-integration 
tricliP^ to shift Pn onto the real axis, giving 


density 

b{Qo,Po,t) = J dP' J dQ' 

X e^“*H(Q) (10) 

where the primes denote integration over all modes ex¬ 
cept Po and Qo- Differentiation with respect to t, appli¬ 
cation of Eq. Q , and integration by parts gives 
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In the usual way of mean-field dynamics ,1^^ we can split 
the force on the centroid into 
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(and we have used the t = 0 contour-integration trick to 
get to the second line), 


cfJiO) = ^JdPjdQ A(Q)P(Q)e-^«"(P'Q) (9) 


Z{Qo) = J dP' J dQ' (15) 


which now contains the (real) ring-polymer 
distribution,E2l and hence no longer suffers from 
the sign problem. Unfortunately, this trick does not 
work for t > 0 (see Sec. IV), so we are stuck with Eq. (|^, 
which motivates us to hnd approximations to Matsubara 
dynamics. 


and Pfluct(Q) is the huctuation force (dehned by Eq. (13) 


as the difference between the exact force and FoiQo))- 
Equation © then splits into 
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III. CENTROID MEAN-FIELD 
APPROXIMATION 

This approximation can be made if A(Q) is a function 
of just the centroid Qo (or Po)P in which case we need 
only the Matsubara dynamics of the centroid reduced 
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This type of expression is encountered in coarse-graining, 
where the integral is sometimes approximated by a gen- 
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IV. ANALYTIC CONTINUATION AT t > 0 


We now return to Eq. Q, which is just Eq. ^ rewrit¬ 
ten in terms of (P,Q) Expressing Cm in terms of these 
coordinates gives 
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FIG. 1. Comparisons of Matsubara, CMD, RPMD and (ex¬ 
act) quantum Kubo-transformed autocorrelation functions, 
for the quartic potential V{q) = q^/4, with mass m = 1, at 
temperature j3 = 2 (in atomic units). The position autocor¬ 
relation functions in (a) are taken from ref. The position- 
squared autocorrelation functions in (b) were calculated nu¬ 
merically using theprocedure described in ref. [2l with M = 7 
Matsubara modes.*^ The differences between the Matsubara 
and exact quantum results show the importance of real-time 
quantum coherence in this model system, the neglect of which 
(in the Matsubara calculations) has blue-shifted and broad¬ 
ened the spectrum. 


eralized Langevin term.^ It is an exact rewriting of 
Eq. 0- Neglect of the integral term gives the mean- 
field approximation 
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is the RPMD Liouvillian (corresponding to the ring- 
polymer hamiltonian i?M(P,Q)) and 



Note that the complete Liouvillian Cm does not corre¬ 
spond to a Hamiltonian in (P, Q) (because the trans¬ 
formation from (P, Q) to (P,Q)is non-canonical), and 
that any resemblance to RPMEpEl is at this stage illu¬ 
sory, since the imaginary parts of n 0, contribute 
terms that cancel the spring terms in . 

If we now try to shift P„, n 0, onto the real axis, 
we find that the dynamics generated by Cm propagates 
an initial distribution of real phase-space points into the 
complex plane, along unstable trajectories. We do not 
know whether the contour-integration trick remains valid 
for such trajectories; even if it does, they appear to be at 
least as difficult to treat numerically as the sign problem 
in Eq. 0. 

However, it is possibl^l^ to follow a path along which 
one gradually moves P„, n y^ 0, towards the real axis 
whilst gradually discarding C^^, such that the dynamics 
remains stable (and the contour-integration trick remains 
valid) at every point along the path. At the end of the 
path, T™ has been completely discarded, and Pn has 
reached the real axis. This results in the approximation, 

C'i^'(t) ~ ^ JdP JdQ 

( 21 ) 


which is CMD.ISIIISI 

Thus CMD corresponds to approximating Matsubara 
dynamics by leaving out the fluctuation term in Eq. (161. 
This result is not a surprise, and is consistent with previ¬ 
ous numerical finding^^that CMD causes errors through 
neglect of fluctuations (see Sec. V). What is new is that 
Eq. (16) gives an explicit formula for these fluctuations, 
in the case that the quantum dynamics can be approxi¬ 
mated by Matsubara dynamics. 


which is RPMD.I^lfiliSl 

A harmonic analysiJSI shows that the main effect of 
discarding is erroneously to shift the Matsubara 
fluctuation frequencies to the ring-polymer frequencies. 
Since C^^ does not act directly on Qq, it follows that 
an RPMD time-correlation function involving linear op¬ 
erators (for which B{Q) = Qo or Pq) will agree initially 
with the Matsubara result, but will then lose accuracy 
as the errors in the fluctuation dynamics couple to the 
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CMD 


RPMD 


satisfies detailed balance, because the centroid mean-field 
force is decoupled from the Matsubara fluctuations 


satisfies detailed balance, because and /ijjj [in 

Eq. (181] independently satisfy detailed balance 


is the centroid mean-field approximation to Matsubara 
dynamics 


has the same centroid mean-field approximation as 
Matsubara dynamics, namely CMD 


is exact for linear TCFs in the harmonic limit, since the 
centroid mean-field force is then equal to the Matsubara 
force 


is exact for linear TCFs in the harmonic limit, since the 
neglected term does not act on the centroid 


gives the exact centroid-averaged Matsubara Liouvillian 
dynamics at t = 0 

suffers from the curvature problem in vibrational spectra 
because of the neglect of the Matsubara fluctuations 

gives the mean-field-averaged Matsubara force on the 
centroid 

breaks down completely for non-linear A and 13 (see 
Fig. lb) because j4(Q) and i3(Q) depend on non-centroid 
modes 


gives the exact Matsubara Liouvillian dynamics at t = 0 

suffers from spurious resonances in vibrational spectra 
because the neglect of shifts the fluctuation frequencies 

gives the exact Matsubara force on the centroid 

breaks down more rapidly for non-linear (than for linear) A 
and B (see Fig. lb) because the neglected term /ijjj acts 
directly on the non-centroid modes 


TABLE 1. Properties of CMD and RPMD derived from Matsubara dynamics (TCF = time-correlation function). 


centroid through the anharmonicity in C/m(Q)- This re¬ 
sult is not a surprise, as the ring-polymer freq uenci es are 
known to interfere with the centroid dynamicsWhat 
is new is that we have identified the approximation made 
by RPMD, namely the neglect of 


V. DISCUSSION 

We have shown that both CMD and RPMD are ap¬ 
proximations to Matsubara dynamics, which, as men¬ 
tioned in the Introduction, is probably the correct way 
to describe quantum statistics and classical dynamics. 
CMD neglects the Matsubara fluctuation term; RPMD 
neglects part of the Matsubara Liouvillian. So far as we 
can tell, there is no direct physical justification that can 
be given for either of these approximations. CMD and 
RPMD are useful because, as has long been known,l2Hsl 
they preserve detailed balance, and satisfy a number of 
important limits. These propertieJi^ (and a few oth¬ 
ers) can be rederived from Matsubara dynamics, and are 
listed in Table I. Note also that CMD and RPMD give 
the same t = 0 leading-order error terms when compared 
with Mat subar a dynamics as with the exact quantum 

dynamics .^^nUI] 

One new finding, less drastic than it first appears, is 
that both CMD and RPMD give qualitatively wrong fluc¬ 
tuation dynamics at barriers. In Matsubara dynamics, 
some of the distributions in ^(t) stretch indefinitely over 
the barrier top, such that a proportion of the distribu¬ 
tion ends up on either side. In CMD and RPMD, all 
of the distribution ends up on one side of the barrier 
(because CMD decouples the fluctuation modes neces¬ 


sary for stretching over the barrier, and RPMD shifts 
the frequencies of these modes from imaginary to reaP^ . 
However, CMD and RPMD are still powerful tools for 
estimating quantum reaction rates, as the exact t = 0 
behaviour of these methods (see Table I) ensures that 
classical rate theory (in the mean-field centroid or ring- 
polymer space) gives lower bound estimates of the 
quantum transition-state theory rate,l^ for the special 
case of a centroid dividing-surface (CMD), and for the 
general case (RPMD). 

The main new result of this work is that, in relating 
CMD and RPMD to Matsubara dynamics, we have ob¬ 
tained explicit formulae for the terms that are left out, 
which may lead to improvements in these methods. For 
example, it might be possible to include approximately 
the Matsubara fluctuation term of Eq. ( [T^ which is miss¬ 
ing in CMD, or to exploit the property that RPMD gives 
the correct Matsubara force on the centroid.!^ 
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APPENDIX: MATSUBARA MODES 

The set of M Matsubara modes Q is defined a^ 

1 ^ 

= lim n = 0, ±1,..., ±(M - l)/2 

N^oo \ N ^' 

^ 1 = 1 

(22) 

where M is odcP and satisfies M <C A; q = {qi},l = 
1,... ,N, are a set of discrete path-integral coordinates 
distributed at equally spaced intervals j3h/N of imaginary 
time, and 

n = 0 

Tir,=ly^^sui{2'Kln/N) n=l,...,(M-l)/2 (23) 

iv'^cos(27rZn/A) n = -l,...,-(M-l)/2 

The momentum coordinates P are similarly defined in 
terms of p. Qq and Pq are the position and momen¬ 
tum centroid coordinates. We define the associated Mat¬ 
subara frequencies = 2mr//3H such that they include 
the sign of n, which gives 0 m(P, Q) the simple form of 
Eq. (|^. 

The functions A(Q) and i3(Q) in Eq. ^ are obtained 
by making the substitutions 

(M-l)/2 

qi = VN (24) 

ra=-(M-l)/2 

into the functions 

TV TV 

^(q)=]^E^(«)’ B{<i) = -Y,B{qi) (25) 

1=1 1=1 

The Matsubara potential Um{Q) is obtained similarly, 
by subsituting for qi in the ring-polymer potential 

1 ^ 

UNiq) = j;^Y.^{qi) (26) 

1=1 

We emphasise that the formulae above and in Sec. II re¬ 
sult from just one approximation, namely decoupling the 
Matsubara modes from the non-Matsubara modes in the 
exact quantum Liouvillian (which causes all Liouvillian 
terms 0{lP) to vanish).!^ 
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I. ANALYTIC CONTINUATION 


To move to the real axis at t = 0, we apply the standard contour-integration trick 
used to take the Fourier tranform of a Gaussian: We construct a rectangular contour with 
top and bottom sides running along ImP„ = —mUnQ-n and the real axis, and vertical 
sections at ReP^ = ±L, then take the limit L —)■ oo. The vertical sections give zero, and 
there are no poles inside the contour, so the integral is unchanged on shifting P„ to the real 
axis. 

We do not know if this trick is valid at t > 0, since the dynamics propagates P and Q 
into the complex plane, resulting in unstable trajectories, which may cause e^’^^B{Q) to 
diverge.^ However, the dynamics does remain stable, and the analytic continuation valid, 
if one gradually discards Pjj]- as one pushes P„ towards the real axis. One can do this by 
letting A change smoothly from 1 to 0 in the Liouvillian 

£, = 4f>+a£!;j (SI) 


whilst setting P„ = H^ — iXmUnQ-n (where n„ and Q-n are real). Writing Cx in terms of 
n„ and Q-m we obtain 


Cx = 


(M-l)/2 
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n=-{M-l)/2 


--i- 
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dUMiQ) 

dQn 


+ m(l - X‘^)ulQn 


d 
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(S2) 


which shows that the dynamics of P and Q maps onto a real dynamics in 11 and Q at every 
value of A between 1 and 0, and thus avoids the unstable trajectories in the complex plane. 
Note that the dynamics generated by Cx satisfies detailed balance (because and 
independently satisfy detailed balance). 
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II. ANALYSIS OF HARMONIC FLUCTUATION TERMS 


If the potential V{q) is harmonic, then Pn in Eq. (7) of the article can be analytically 
continued to the real axis for f > 0, giving the Matsubara equations of motion 

Qn(t) = Qn cos ut H- — siu ult + i—Q_n sin uit 

mu u 

_2 

Pn(t) = Pn COS ut - siu Wt — i — P_„ siu (S3) 

u u 

where Un = + u‘^ are the ring-polymer frequencies. On discarding we obtain the 

RPMD equations of motion 

~ ~ — T^n — 

Qn{t) = Qn COS Unt H-^ siu Unt 

mUn 

Pn(t) = Pn sin Unt — mUnQn^^^^nt (S4) 

Hence the omission of changes the fluctuation frequencies u to the ring-polymer frequen¬ 
cies Un- Curiously, the amplitude of the fluctuations is unchanged, since for both Matsubara 
and RPMD dynamics, the (M — l)/2 quantities, 

(Qu + Q-n^ + -Pn + ^-n (S5) 

are constants of the motion. This analysis can be repeated for a parabolic barrier, where it 
is found that discarding changes the fluctuation frequencies from the imaginary barrier 
frequency i\u\ to the real ring-polymer frequencies Un = — \u\^ (if jSh < 2 t^/\u\). 
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